## This script produces Goodness of Fit diagnostics for "The Coevolution of Networks of Interstate Support, Interstate Threat and Civil War" 
## Since the estimation involves simulations, results may vary slightly from the published results.


library(reshape2)
library(network)
library(sna)
library(RSiena)
library(ggplot2)
library(gdata)
library(countrycode)
library(parallel)
library(ROCR)
library(PRROC)
library(miceadds)
library(sandwich)

rm(list=ls(all=TRUE))

#dir <- "[Directory Information Here]"
dir <- "C:/Users/rapiduser/Box/Relation creation/JOPreplication/"

setwd(dir)

(n.clus <- detectCores() - 1)


###############################################################################
##OOS GOF
###############################################################################

## The starting point for this script came from Kinne & Bunte, "Guns or Money," BJPS 2018

keep(dir, n.clus, sure=TRUE)

nets.total <- read.csv("relation_network.csv", header=TRUE, row.names=1)
mons.total <- read.csv("relation_nodal.csv", header=TRUE, row.names=1)

## Set some parameters.
## We're using data only through the trainyr, and saving testyr for prediction

#############################################################################################################
#2004 training period
##############################################################################################################
trainyr<-2004

yr <- seq(1980, trainyr)
yr.lo <- yr[1:(length(yr)-1)]
compchanges <- read.csv("compchanges04.csv", header=TRUE, row.names=1)

#Remove country 341, which enters in 2006
nets.total <- nets.total[nets.total$ccode1!=341 & nets.total$ccode2!=341,]
mons.total <- mons.total[mons.total$ccode1!=341,]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

## Make composition change list 
write.table(compchanges, "cc.dat",row.names=FALSE, col.names=FALSE,
            quote=FALSE)
changes <- sienaCompositionChangeFromFile("cc.dat")
names(changes) <- id

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not substracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Changing dyadic covariates, omitting the constant ones
dyad.effects <- d.cons[ !d.cons %in% c("contiguous") ]
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), varDyadCovar(array(do.call("c", get(x)), dim=c(N,N,nlo))))
}; rm(x)

## Constant dyadic covariates
for (x in d.cons[ d.cons %in% c("contiguous")]) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), varCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net,
  
  changes
  
)

eff <- getEffects(netdata)

###Start with base model, and then add in network effects

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)

## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                       interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                       interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect
eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                       interaction1="support.net", include=T, type="eval")

##Now add information from the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")

## Estimate models and save the results as the training set
setwd(paste(dir, "/Output", sep=""))

tm <- Sys.time()
model <- sienaAlgorithmCreate(useStdInits=F, projname="Estimates.GOF.train.04",
                              nsub=4, n3=5000, maxlike=F, firstg=.03,
                              modelType=c(support.net=1,threat.net=1))
output <- siena07(model,data=netdata,effects=eff,batch=TRUE,
                  verbose=FALSE,useCluster=TRUE,nbrNodes=n.clus,returnDeps=FALSE)
save(output, file="output.GOF.train.paper_04")

rm(output)




######################################################################################
## Now specify models for the OOS periods
######################################################################################

###############################
##2005
###############################

keep(dir, n.clus, mons.total, nets.total, trainyr, sure=TRUE)
setwd(dir)

testyr<-2005

## Specify data for out-of-sample prediction
yr <- seq((testyr-1), testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Must specify all covariates as constant (i.e., not changing over time)
dyad.effects <- d.cons
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}; rm(x)

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), coCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
  
)

eff <- getEffects(netdata)

###Specify the SAOM model

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect
eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add characteristics of the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")


## Grab estimates from training sample
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.train.paper_04")

## Save rate and degree start values, since these must differ
## from training period.
rp1 <- eff[eff$effectName=="basic rate parameter support.net",]$initialValue
rp2 <- eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp3 <- eff[eff$effectName=="basic rate parameter threat.net",]$initialValue
rp4 <- eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp5 <- eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue
rp6 <- eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue

## Update thetas with estimates from training period
eff <- updateTheta(eff,output)
## Then replace rate and density parameters
eff[eff$effectName=="basic rate parameter support.net",]$initialValue <- rp1
eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp2
eff[eff$effectName=="basic rate parameter threat.net",]$initialValue <- rp3
eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp4
eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue <- rp5
eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue <- rp6

## No subroutines in stage 2, since we're not estimating parameters. Instead,
## just simulate the network a few thousand times.
model.SIM <- sienaModelCreate(useStdInits = F,
                              projname="Simulation-GOF-Predict-05",
                              nsub=0, n3=3000, maxlike=F, cond=F, modelType=c(support.net=1,threat.net=1))

## Call the estimation/simulation function. Be sure to return the dependent
## simulated matrices.
SIM.out.predict <- siena07(model.SIM, data=netdata, effects=eff,
                           batch=TRUE, verbose=FALSE, useCluster=TRUE,
                           nbrNodes=n.clus, initC=TRUE, returnDeps=TRUE)
save(SIM.out.predict,file=paste("output.GOF.predict.paper_05"))


keep(mons.total, nets.total, dir, n.clus, id, N, yr, Ss.NET, Ts.NET, intra_intens, trainyr, testyr, sure=TRUE)

## Load SAOM predictions
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.predict.paper_05")



## Collect SAOM predictions for threat and intrastate conflict
sims <- SIM.out.predict$sims

#Start for threat ties
## Grab predicted SAOM networks for last time period
base <- matrix(0,N,N,dimnames=list(id,id))
for ( i in 1:length(sims) ) {
        adj <- matrix(0,N,N)
        ## The final year is the one we're predicting
        edges <- sims[[i]][[1]][["threat.net"]][[(length(yr)-1)]]
        adj[edges[, 1:2]] <- edges[,3]
        base <- base + adj
        rm(adj,edges)
}
## Weight total predicted ties by number of simulated networks
base <- base/length(sims); diag(base) <- NA

## Get the * observed* networks for last time period
obs <- Ts.NET[[length(yr)]]
diag(obs) <- NA
    
## Combine predicted and observed vectors
saom.pre.threat.net.05 <- na.omit(cbind(as.data.frame(as.table(base)),as.data.frame(as.table(obs))[,3]))
names(saom.pre.threat.net.05) <- c("ccode1","ccode2","predict","observe")    

#Now for Intrastate violence -- specifically looking at all armed conflict
base<-as.vector(matrix(0,nrow=N,dimnames=list(id)))
for ( i in 1:length(sims) ) {
  ## The final year is the one we're predicting
  ac <- sims[[i]][[1]][["intra_intens.net"]][[(length(yr)-1)]]
  ac[ac>=1]<-1
  base <- base + ac
  rm(ac)
}

## Weight total predicted ties by number of simulated networks
base <- base/length(sims)

## Get the observed civil wars from last time period
intra.obs<-intra_intens[,length(yr)]

##combine predicted and observed
op <- na.omit(cbind(base,intra.obs,id))
saom.pre.intra_intens.05<-as.data.frame(op)
names(saom.pre.intra_intens.05) <- c("predict","observe","ccode1")
saom.pre.intra_intens.05$observe[saom.pre.intra_intens.05$observe>=1]<-1

rm(base,op,sims)


#################################################

## Now estimate a logit model, predict out of sample, and compare
## the logit to the SAOM predictions

## Reorganize the dataset for the logit model. Use the period up until the test year.
## We're going to grab all the same data as for the SAOMs, then
## restructure it into a data frame.
yr <- seq(1980, testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with the DVs
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))

#########################################################
## Then do all the covariates for both networks
#########################################################


## Names of dyadic and monadic covariates, plus behavior variable
#Create new name for the intra_intens that will be a covariate
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]
mons.total$intra_intensc<-mons.total$intra_intens
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year") ]


## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################################
## Now convert all the data matrices to a data frame - one dyadic and one monadic
####################################################

dat <- lapply(Ss.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
    dat[[w]]$year <- yr[w]
}
dat <- do.call("rbind", dat)
names(dat) <- c("ccode1","ccode2","Ss","year")

tmp <- lapply(Ts.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
    tmp[[w]]$year <- yr[w]
}
tmp <- do.call("rbind", tmp)
names(tmp) <- c("ccode1","ccode2","Ts","year")

dat <- merge(dat,tmp); rm(tmp)

## Add all the covariates to the data frame, dyadic covariates first. 
##Lag these all by a year by treating them as the t+1 data and then merge with the DVs.
for (w in d.cons) {
    tmp <- lapply(get(w), function(x) as.data.frame(as.table(x)))
    for (z in 1:length(tmp)) {
        tmp[[z]]$year <- yr[z] + 1
    }
    tmp <- do.call("rbind", tmp)
    names(tmp) <- c("ccode1","ccode2",w,"year")
    dat <- merge(dat,tmp,all.x=TRUE)
    rm(tmp)
}; rm(w)

## Now merge monadic covariates.
for (w in m.cons) {
    tmp <- as.data.frame(as.table(get(w)))
    names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
    tmp$year <- as.numeric(as.character(tmp$year)) + 1
    dat <- merge(dat,tmp,all.x=TRUE)
    names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
    dat <- merge(dat,tmp,all.x=TRUE)
    if ( w == "dem" ) {
      dat[[paste(w,".both",sep="")]] <- dat[[paste(w,".cc1", sep="")]] * dat[[paste(w,".cc2", sep="")]]
   }
    rm(tmp)
}

for (w in b.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}

## Drop "self dyads
dat <- dat[dat$ccode1 != dat$ccode2,]

####Make the monadic data for models of intrastate conflict
dat_mon <-mons.total[,c("ccode1","year","intra_intens")]
names(dat_mon) <- c("ccode1","year","intra_intens")
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",w)
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat_mon <- merge(dat_mon,tmp,all.x=TRUE)
  rm(tmp)
}


## Make unique identifiers
dat$dyadid <- (as.numeric(dat$ccode1) * 1000) + as.numeric(dat$ccode2)

## Now collect all the nonmissing observations
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
## Specify training set and validation set
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)

##### Estimate logit models 
logit.threatn.train <- glm(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.threatn.train, dt[idx.out,], type="response", )
pred.use.threatn.05 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.threatn.05) <- c("observe","ccode1","ccode2","predict")

logit.basic.threatn.train <- glm(Ts ~ intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.basic.threatn.train, dt[idx.out,], type="response", )
pred.use.basic.threatn.05 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.basic.threatn.05) <- c("observe","ccode1","ccode2","predict")



## Now do the same for Intrastate conflict, constraining intra_intens to binary
## Specify training set and validation set
dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)

logit.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.intran.05 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.intran.05) <- c("observe","ccode1","predict")

logit.basic.intran.train <- glm(intra_intens ~  intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.basic.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.basic.intran.05 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.basic.intran.05) <- c("observe","ccode1","predict")

#####Generate coef plots of the logit models, clustering the SEs, Appendix Figure 11
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)
logit.threatn.train <-glm.cluster(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, cluster='dyadid', data=dt, family=binomial(link="logit"), subset=idx.in)

dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)
logit.intran.train <- glm.cluster(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, cluster='ccode1', family=binomial(link="logit"), subset=idx.in)


## First do threats logit
output.threat <- data.frame(summary(logit.threatn.train))
out.threat <- cbind(output.threat$Estimate,output.threat$Std..Error)
rownames(out.threat) <- rownames(output.threat)
for (i in 1:nrow(out.threat)) { # Rescale means and SEs for forest plot with all SEs at 0.5
  ad <- 0.5/out.threat[i,2]
  out.threat[i,1] <- out.threat[i,1] * ad
  out.threat[i,2] <- out.threat[i,2] * ad
  rm(ad)
}; rm(i)
#Rescale intercept and reciprocity
out.threat[1,1] <- out.threat[1,1] * 0.1
out.threat[1,2] <- out.threat[1,2] * 0.1
out.threat[17,1] <- out.threat[17,1] * 0.1
out.threat[17,2] <- out.threat[17,2] * 0.1
threat.names <-  c("Intercept (x 0.1)","Support Cross Reciprocity",
                   "i's Intrastate Conflict Intensity","j's Intrastate Conflict Intensity",
                   "Common IGO Memberships","Ideal Point Distance","Contiguity",
                   "i's Democracy","j's Democracy","Joint Democracy",
                   "i's Military Expenditures Ratio","j's Military Expenditures Ratio",
                   "i's Energy Consumption Per Capita","j's Energy Consumption Per Capita",
                   "i's Military Capabilities Change","j's Military Capabilities Change",
                   "Reciprocity (x 0.1)",
                   "i's Support Military Spending","i's Support Ideal-Point Closeness",
                   "j's Support Military Spending","j's Support Ideal-Point Closeness",
                   "i's Support IP Closeness x Military Spending","j's Support IP Closeness x Military Spending")
dthreat <- data.frame(
  nms=threat.names,
  x=out.threat[,1],
  xlo=(out.threat[,1] - (1.96 * out.threat[,2])),
  xhi=(out.threat[,1] + (1.96 * out.threat[,2]))
)
p.cols <- vector()
for (ii in 1:nrow(dthreat)) {
  if (dthreat[ii,"xlo"] < 0 && dthreat[ii,"xhi"] > 0) {
    p.cols[ii] <- "#00b2a6" # Color for insignificant estimates
  } else {
    p.cols[ii] <- "#8519d7" # Color for significant estimates
  }
}; rm(ii)
dthreat$p.cols <- p.cols
#Adjust the order
dthreat <- dthreat[ c(1,17,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,22,19,20,23,21), ]

marg <- c(0.1,0.5,0.1,0)

p2 <- ggplot(dthreat) + geom_hline(yintercept=0, lty=2) +
  geom_linerange(mapping=aes(x=seq(1,nrow(dthreat)), ymin=xlo, ymax=xhi), color=dthreat$p.cols, linetype=1, size=1, alpha=0.5) +
  geom_point(mapping=aes(y=x, x=seq(1,nrow(dthreat))), size=3, shape=16, color=dthreat$p.cols) +
  scale_x_continuous(breaks=seq(1:nrow(dthreat)), labels=as.character(dthreat$nms)) +
  xlab("") + ylab("Rescaled estimates + 95% CIs") + ggtitle("Logit Model of Threat") + coord_flip() + 
  theme(
    panel.grid.minor.y=element_blank(),
    axis.text.y=element_text(color="black"),
    plot.margin=unit(marg,"cm"),
    plot.title = element_text(hjust = 0.5)
  )

setwd(paste(dir, "/Figures", sep=""))
pdf("ThreatEstimates_paper_logit.pdf",width=6,height=8,paper="special")
par(mar=c(0,0,0,0))
plot(p2)
dev.off()


#Now Intrastate Logit
output.intra <- data.frame(summary(logit.intran.train))
out.intra <- cbind(output.intra$Estimate,output.intra$Std..Error)
rownames(out.intra) <- rownames(output.intra)
for (i in 1:nrow(out.intra)) { # Rescale means and SEs for forest plot with all SEs at 0.5
  ad <- 0.5/out.intra[i,2]
  out.intra[i,1] <- out.intra[i,1] * ad
  out.intra[i,2] <- out.intra[i,2] * ad
  rm(ad)
}; rm(i)
intra.names <-  c("Intercept","i's Intrastate Conflict Intensity",
                  "i's Intrastate Conflict Intensity (sq)","i's Democracy",
                  "i's Military Expenditures Ratio","i's Energy Consumption Per Capita",
                  "i's Military Capabilities Change","i's Threat Military Spending",
                  "i's Support Military Spending",
                  "i's Support Ideal-Point Closeness",
                  "i's Support IP Closeness x Military Spending")
dintra <- data.frame(
  nms=intra.names,
  x=out.intra[,1],
  xlo=(out.intra[,1] - (1.96 * out.intra[,2])),
  xhi=(out.intra[,1] + (1.96 * out.intra[,2]))
)
p.cols <- vector()
for (ii in 1:nrow(dintra)) {
  if (dintra[ii,"xlo"] < 0 && dintra[ii,"xhi"] > 0) {
    p.cols[ii] <- "#00b2a6" # Color for insignificant estimates
  } else {
    p.cols[ii] <- "#8519d7" # Color for significant estimates
  }
}; rm(ii)
dintra$p.cols <- p.cols
#Adjust the order
dintra <- dintra[ c(1,2,3,4,5,6,7,8,9,11,10), ]

marg <- c(0.1,0.5,0.1,0)

p2 <- ggplot(dintra) + geom_hline(yintercept=0, lty=2) +
  geom_linerange(mapping=aes(x=seq(1,nrow(dintra)), ymin=xlo, ymax=xhi), color=dintra$p.cols, linetype=1, size=1, alpha=0.5) +
  geom_point(mapping=aes(y=x, x=seq(1,nrow(dintra))), size=3, shape=16, color=dintra$p.cols) +
  scale_x_continuous(breaks=seq(1:nrow(dintra)), labels=as.character(dintra$nms)) +
  xlab("") + ylab("Rescaled estimates + 95% CIs") + ggtitle("Logit Model of Intrastate Conflict") + coord_flip() + 
  theme(
    panel.grid.minor.y=element_blank(),
    axis.text.y=element_text(color="black"),
    plot.margin=unit(marg,"cm"),
    plot.title = element_text(hjust = 0.5)
  )

setwd(paste(dir, "/Figures", sep=""))
pdf("Intra_IntensEstimates_paper_logit.pdf",width=6,height=4,paper="special")
par(mar=c(0,0,0,0))
plot(p2)
dev.off()



###############################
##2006
###############################

keep(dir, n.clus, mons.total, nets.total, trainyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, sure=TRUE)
setwd(dir)

testyr<-2006

## Specify data for out-of-sample prediction
yr <- seq((testyr-1), testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with making the DV network matrices
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Must specify all covariates as constant (i.e., not changing over time)
dyad.effects <- d.cons
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}; rm(x)

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), coCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
  
)

eff <- getEffects(netdata)

###Specify the SAOM model

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect

eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add information from the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")


## Grab estimates from training sample
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.train.paper_04")

## Save rate and degree start values, since these must differ
## from training period.
rp1 <- eff[eff$effectName=="basic rate parameter support.net",]$initialValue
rp2 <- eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp3 <- eff[eff$effectName=="basic rate parameter threat.net",]$initialValue
rp4 <- eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp5 <- eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue
rp6 <- eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue

## Update thetas with estimates from training period
eff <- updateTheta(eff,output)
## Then replace rate and density parameters
eff[eff$effectName=="basic rate parameter support.net",]$initialValue <- rp1
eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp2
eff[eff$effectName=="basic rate parameter threat.net",]$initialValue <- rp3
eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp4
eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue <- rp5
eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue <- rp6

## No subroutines in stage 2, since we're not estimating parameters. Instead,
## just simulate the network a few thousand times.
model.SIM <- sienaModelCreate(useStdInits = F,
                              projname="Simulation-GOF-Predict-06",
                              nsub=0, n3=3000, maxlike=F, cond=F, modelType=c(support.net=1,threat.net=1))

## Call the estimation/simulation function. Be sure to return the dependent
## simulated matrices.
(n.clus <- detectCores() - 1)
SIM.out.predict <- siena07(model.SIM, data=netdata, effects=eff,
                           batch=TRUE, verbose=FALSE, useCluster=TRUE,
                           nbrNodes=n.clus, initC=TRUE, returnDeps=TRUE)
save(SIM.out.predict,file=paste("output.GOF.predict.paper_06"))


keep(mons.total, nets.total, dir, n.clus, id, N, yr, Ss.NET, Ts.NET, intra_intens, trainyr, testyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, sure=TRUE)

## Load SAOM predictions
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.predict.paper_06")



## Collect SAOM predictions for threat and intrastate conflict
sims <- SIM.out.predict$sims

#Start for threat ties
## Grab predicted SAOM networks for last time period
base <- matrix(0,N,N,dimnames=list(id,id))
for ( i in 1:length(sims) ) {
  adj <- matrix(0,N,N)
  ## The final year is the one we're predicting
  edges <- sims[[i]][[1]][["threat.net"]][[(length(yr)-1)]]
  adj[edges[, 1:2]] <- edges[,3]
  base <- base + adj
  rm(adj,edges)
}
## Weight total predicted ties by number of simulated networks
base <- base/length(sims); diag(base) <- NA

## Get the * observed* networks for last time period
obs <- Ts.NET[[length(yr)]]
diag(obs) <- NA

## Combine predicted and observed vectors
saom.pre.threat.net.06 <- na.omit(cbind(as.data.frame(as.table(base)),as.data.frame(as.table(obs))[,3]))
names(saom.pre.threat.net.06) <- c("ccode1","ccode2","predict","observe")    

#Now for Intrastate violence -- specifically looking at all armed conflict
base<-as.vector(matrix(0,nrow=N,dimnames=list(id)))
for ( i in 1:length(sims) ) {
  ## The final year is the one we're predicting
  ac <- sims[[i]][[1]][["intra_intens.net"]][[(length(yr)-1)]]
  ac[ac>=1]<-1
  base <- base + ac
  rm(ac)
}

## Weight total predicted ties by number of simulated networks
base <- base/length(sims)

## Get the observed civil wars from last time period
intra.obs<-intra_intens[,length(yr)]

##combine predicted and observed
op <- na.omit(cbind(base,intra.obs,id))
saom.pre.intra_intens.06<-as.data.frame(op)
names(saom.pre.intra_intens.06) <- c("predict","observe","ccode1")
saom.pre.intra_intens.06$observe[saom.pre.intra_intens.06$observe>=1]<-1

rm(base,op,sims)


#################################################

## Now estimate a logit model, predict out of sample, and compare
## the logit to the SAOM predictions

## Reorganize the dataset for the logit model. Use the period up until the test year.
## We're going to grab all the same data as for the SAOMs, then
## restructure it into a data frame.
yr <- seq(1980, testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with the DVs
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))

#########################################################
## Then do all the covariates for both networks
#########################################################


## Names of dyadic and monadic covariates, plus behavior variable
#Create new name for the intra_intens that will be a covariate
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]
mons.total$intra_intensc<-mons.total$intra_intens
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year") ]


## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################################
## Now convert all the data matrices to a data frame - one dyadic and one monadic
####################################################

dat <- lapply(Ss.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  dat[[w]]$year <- yr[w]
}
dat <- do.call("rbind", dat)
names(dat) <- c("ccode1","ccode2","Ss","year")

tmp <- lapply(Ts.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  tmp[[w]]$year <- yr[w]
}
tmp <- do.call("rbind", tmp)
names(tmp) <- c("ccode1","ccode2","Ts","year")

dat <- merge(dat,tmp); rm(tmp)

## Add all the covariates to the data frame,
## dyadic covariates first. Lag these all
## by a year.
for (w in d.cons) {
  tmp <- lapply(get(w), function(x) as.data.frame(as.table(x)))
  for (z in 1:length(tmp)) {
    tmp[[z]]$year <- yr[z] + 1
  }
  tmp <- do.call("rbind", tmp)
  names(tmp) <- c("ccode1","ccode2",w,"year")
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}; rm(w)

## Now merge monadic covariates.
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  if ( w == "dem" ) {
    dat[[paste(w,".both",sep="")]] <- dat[[paste(w,".cc1", sep="")]] * dat[[paste(w,".cc2", sep="")]]
  }
  rm(tmp)
}

for (w in b.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}

## Drop "self dyads
dat <- dat[dat$ccode1 != dat$ccode2,]

####Make the monadic data for models of intrastate conflict
dat_mon <-mons.total[,c("ccode1","year","intra_intens")]
names(dat_mon) <- c("ccode1","year","intra_intens")
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",w)
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat_mon <- merge(dat_mon,tmp,all.x=TRUE)
  rm(tmp)
}


## Make unique identifiers
dat$dyadid <- (as.numeric(dat$ccode1) * 1000) + as.numeric(dat$ccode2)

## Now collect all the nonmissing observations
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
## Specify training set and validation set
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)



##### Estimate logit models
logit.threatn.train <- glm(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.threatn.train, dt[idx.out,], type="response", )
pred.use.threatn.06 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.threatn.06) <- c("observe","ccode1","ccode2","predict")

logit.basic.threatn.train <- glm(Ts ~ intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.basic.threatn.train, dt[idx.out,], type="response", )
pred.use.basic.threatn.06 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.basic.threatn.06) <- c("observe","ccode1","ccode2","predict")



## Now do the same for Intrastate conflict, constraining intra_intens to binary
## Specify training set and validation set
dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)

logit.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.intran.06 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.intran.06) <- c("observe","ccode1","predict")

logit.basic.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.basic.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.basic.intran.06 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.basic.intran.06) <- c("observe","ccode1","predict")


###############################
##2007
###############################

keep(dir, n.clus, mons.total, nets.total, trainyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, sure=TRUE)
setwd(dir)

testyr<-2007

## Specify data for out-of-sample prediction
yr <- seq((testyr-1), testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Must specify all covariates as constant (i.e., not changing over time)
dyad.effects <- d.cons
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}; rm(x)

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), coCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
  
)

eff <- getEffects(netdata)

###Specify the SAOM model

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)

## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect

eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add information from the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")


## Grab estimates from training sample
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.train.paper_04")

## Save rate and degree start values, since these must differ
## from training period.
rp1 <- eff[eff$effectName=="basic rate parameter support.net",]$initialValue
rp2 <- eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp3 <- eff[eff$effectName=="basic rate parameter threat.net",]$initialValue
rp4 <- eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp5 <- eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue
rp6 <- eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue

## Update thetas with estimates from training period
eff <- updateTheta(eff,output)
## Then replace rate and density parameters
eff[eff$effectName=="basic rate parameter support.net",]$initialValue <- rp1
eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp2
eff[eff$effectName=="basic rate parameter threat.net",]$initialValue <- rp3
eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp4
eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue <- rp5
eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue <- rp6

## No subroutines in stage 2, since we're not estimating parameters. Instead,
## just simulate the network a few thousand times.
model.SIM <- sienaModelCreate(useStdInits = F,
                              projname="Simulation-GOF-Predict-07",
                              nsub=0, n3=3000, maxlike=F, cond=F, modelType=c(support.net=1,threat.net=1))

## Call the estimation/simulation function. Be sure to return the dependent
## simulated matrices.
(n.clus <- detectCores() - 1)
SIM.out.predict <- siena07(model.SIM, data=netdata, effects=eff,
                           batch=TRUE, verbose=FALSE, useCluster=TRUE,
                           nbrNodes=n.clus, initC=TRUE, returnDeps=TRUE)
save(SIM.out.predict,file=paste("output.GOF.predict.paper_07"))


keep(mons.total, nets.total, dir, n.clus, id, N, yr, Ss.NET, Ts.NET, intra_intens, trainyr, testyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, sure=TRUE)

## Load SAOM predictions
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.predict.paper_07")



## Collect SAOM predictions for threat and intrastate conflict
sims <- SIM.out.predict$sims

#Start for threat ties
## Grab predicted SAOM networks for last time period
base <- matrix(0,N,N,dimnames=list(id,id))
for ( i in 1:length(sims) ) {
  adj <- matrix(0,N,N)
  ## The final year is the one we're predicting
  edges <- sims[[i]][[1]][["threat.net"]][[(length(yr)-1)]]
  adj[edges[, 1:2]] <- edges[,3]
  base <- base + adj
  rm(adj,edges)
}
## Weight total predicted ties by number of simulated networks
base <- base/length(sims); diag(base) <- NA

## Get the * observed* networks for last time period
obs <- Ts.NET[[length(yr)]]
diag(obs) <- NA

## Combine predicted and observed vectors
saom.pre.threat.net.07 <- na.omit(cbind(as.data.frame(as.table(base)),as.data.frame(as.table(obs))[,3]))
names(saom.pre.threat.net.07) <- c("ccode1","ccode2","predict","observe")    

#Now for Intrastate violence -- specifically looking at all armed conflict
base<-as.vector(matrix(0,nrow=N,dimnames=list(id)))
for ( i in 1:length(sims) ) {
  ## The final year is the one we're predicting
  ac <- sims[[i]][[1]][["intra_intens.net"]][[(length(yr)-1)]]
  ac[ac>=1]<-1
  base <- base + ac
  rm(ac)
}

## Weight total predicted ties by number of simulated networks
base <- base/length(sims)

## Get the observed civil wars from last time period
intra.obs<-intra_intens[,length(yr)]

##combine predicted and observed
op <- na.omit(cbind(base,intra.obs,id))
saom.pre.intra_intens.07<-as.data.frame(op)
names(saom.pre.intra_intens.07) <- c("predict","observe","ccode1")
saom.pre.intra_intens.07$observe[saom.pre.intra_intens.07$observe>=1]<-1

rm(base,op,sims)


#################################################

## Now estimate a logit model, predict out of sample, and compare
## the logit to the SAOM predictions

## Reorganize the dataset for the logit model. Use the period up until the test year.
## We're going to grab all the same data as for the SAOMs, then
## restructure it into a data frame.
yr <- seq(1980, testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with the DVs
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))

#########################################################
## Then do all the covariates for both networks
#########################################################


## Names of dyadic and monadic covariates, plus behavior variable
#Create new name for the intra_intens that will be a covariate
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]
mons.total$intra_intensc<-mons.total$intra_intens
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year") ]


## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################################
## Now convert all the data matrices to a data frame - one dyadic and one monadic
####################################################

dat <- lapply(Ss.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  dat[[w]]$year <- yr[w]
}
dat <- do.call("rbind", dat)
names(dat) <- c("ccode1","ccode2","Ss","year")

tmp <- lapply(Ts.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  tmp[[w]]$year <- yr[w]
}
tmp <- do.call("rbind", tmp)
names(tmp) <- c("ccode1","ccode2","Ts","year")

dat <- merge(dat,tmp); rm(tmp)

## Add all the covariates to the data frame,
## dyadic covariates first. Lag these all
## by a year.
for (w in d.cons) {
  tmp <- lapply(get(w), function(x) as.data.frame(as.table(x)))
  for (z in 1:length(tmp)) {
    tmp[[z]]$year <- yr[z] + 1
  }
  tmp <- do.call("rbind", tmp)
  names(tmp) <- c("ccode1","ccode2",w,"year")
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}; rm(w)

## Now merge monadic covariates.
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  if ( w == "dem" ) {
    dat[[paste(w,".both",sep="")]] <- dat[[paste(w,".cc1", sep="")]] * dat[[paste(w,".cc2", sep="")]]
  }
  rm(tmp)
}

for (w in b.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}

## Drop "self dyads
dat <- dat[dat$ccode1 != dat$ccode2,]

####Make the monadic data for models of intrastate conflict
dat_mon <-mons.total[,c("ccode1","year","intra_intens")]
names(dat_mon) <- c("ccode1","year","intra_intens")
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",w)
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat_mon <- merge(dat_mon,tmp,all.x=TRUE)
  rm(tmp)
}


## Make unique identifiers
dat$dyadid <- (as.numeric(dat$ccode1) * 1000) + as.numeric(dat$ccode2)

## Now collect all the nonmissing observations
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
## Specify training set and validation set
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)



##### Estimate logit models
logit.threatn.train <- glm(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.threatn.train, dt[idx.out,], type="response", )
pred.use.threatn.07 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.threatn.07) <- c("observe","ccode1","ccode2","predict")

logit.basic.threatn.train <- glm(Ts ~ intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.basic.threatn.train, dt[idx.out,], type="response", )
pred.use.basic.threatn.07 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.basic.threatn.07) <- c("observe","ccode1","ccode2","predict")

## Now do the same for Intrastate conflict, constraining intra_intens to binary
## Specify training set and validation set
dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)

logit.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.intran.07 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.intran.07) <- c("observe","ccode1","predict")

logit.basic.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.basic.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.basic.intran.07 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.basic.intran.07) <- c("observe","ccode1","predict")


###############################
##2008
###############################

keep(dir, n.clus, mons.total, nets.total, trainyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, saom.pre.threat.net.07, saom.pre.intra_intens.07, pred.use.threatn.07, pred.use.intran.07, pred.use.basic.threatn.07, pred.use.basic.intran.07,sure=TRUE)
setwd(dir)

testyr<-2008

## Specify data for out-of-sample prediction
yr <- seq((testyr-1), testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Must specify all covariates as constant (i.e., not changing over time)
dyad.effects <- d.cons
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}; rm(x)

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), coCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
  
)

eff <- getEffects(netdata)

###Specify the SAOM model

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect

eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add information from the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")


## Grab estimates from training sample
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.train.paper_04")

## Save rate and degree start values, since these must differ
## from training period.
rp1 <- eff[eff$effectName=="basic rate parameter support.net",]$initialValue
rp2 <- eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp3 <- eff[eff$effectName=="basic rate parameter threat.net",]$initialValue
rp4 <- eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp5 <- eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue
rp6 <- eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue

## Update thetas with estimates from training period
eff <- updateTheta(eff,output)
## Then replace rate and density parameters
eff[eff$effectName=="basic rate parameter support.net",]$initialValue <- rp1
eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp2
eff[eff$effectName=="basic rate parameter threat.net",]$initialValue <- rp3
eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp4
eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue <- rp5
eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue <- rp6

## No subroutines in stage 2, since we're not estimating parameters. Instead,
## just simulate the network a few thousand times.
model.SIM <- sienaModelCreate(useStdInits = F,
                              projname="Simulation-GOF-Predict-08",
                              nsub=0, n3=3000, maxlike=F, cond=F, modelType=c(support.net=1,threat.net=1))

## Call the estimation/simulation function. Be sure to return the dependent
## simulated matrices.
(n.clus <- detectCores() - 1)
SIM.out.predict <- siena07(model.SIM, data=netdata, effects=eff,
                           batch=TRUE, verbose=FALSE, useCluster=TRUE,
                           nbrNodes=n.clus, initC=TRUE, returnDeps=TRUE)
save(SIM.out.predict,file=paste("output.GOF.predict.paper_08"))


keep(mons.total, nets.total, dir, n.clus, id, N, yr, Ss.NET, Ts.NET, intra_intens, trainyr, testyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, saom.pre.threat.net.07, saom.pre.intra_intens.07, pred.use.threatn.07, pred.use.intran.07, pred.use.basic.threatn.07, pred.use.basic.intran.07, sure=TRUE)

## Load SAOM predictions
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.predict.paper_08")



## Collect SAOM predictions for threat and intrastate conflict
sims <- SIM.out.predict$sims

#Start for threat ties
## Grab predicted SAOM networks for last time period
base <- matrix(0,N,N,dimnames=list(id,id))
for ( i in 1:length(sims) ) {
  adj <- matrix(0,N,N)
  ## The final year is the one we're predicting
  edges <- sims[[i]][[1]][["threat.net"]][[(length(yr)-1)]]
  adj[edges[, 1:2]] <- edges[,3]
  base <- base + adj
  rm(adj,edges)
}
## Weight total predicted ties by number of simulated networks
base <- base/length(sims); diag(base) <- NA

## Get the * observed* networks for last time period
obs <- Ts.NET[[length(yr)]]
diag(obs) <- NA

## Combine predicted and observed vectors
saom.pre.threat.net.08 <- na.omit(cbind(as.data.frame(as.table(base)),as.data.frame(as.table(obs))[,3]))
names(saom.pre.threat.net.08) <- c("ccode1","ccode2","predict","observe")    

#Now for Intrastate violence -- specifically looking at all armed conflict
base<-as.vector(matrix(0,nrow=N,dimnames=list(id)))
for ( i in 1:length(sims) ) {
  ## The final year is the one we're predicting
  ac <- sims[[i]][[1]][["intra_intens.net"]][[(length(yr)-1)]]
  ac[ac>=1]<-1
  base <- base + ac
  rm(ac)
}

## Weight total predicted ties by number of simulated networks
base <- base/length(sims)

## Get the observed civil wars from last time period
intra.obs<-intra_intens[,length(yr)]

##combine predicted and observed
op <- na.omit(cbind(base,intra.obs,id))
saom.pre.intra_intens.08<-as.data.frame(op)
names(saom.pre.intra_intens.08) <- c("predict","observe","ccode1")
saom.pre.intra_intens.08$observe[saom.pre.intra_intens.08$observe>=1]<-1

rm(base,op,sims)


#################################################

## Now estimate a logit model, predict out of sample, and compare
## the logit to the SAOM predictions

## Reorganize the dataset for the logit model. Use the period up until the test year.
## We're going to grab all the same data as for the SAOMs, then
## restructure it into a data frame.
yr <- seq(1980, testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with the DVs
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))

#########################################################
## Then do all the covariates for both networks
#########################################################


## Names of dyadic and monadic covariates, plus behavior variable
#Create new name for the intra_intens that will be a covariate
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]
mons.total$intra_intensc<-mons.total$intra_intens
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year") ]


## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################################
## Now convert all the data matrices to a data frame - one dyadic and one monadic
####################################################

dat <- lapply(Ss.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  dat[[w]]$year <- yr[w]
}
dat <- do.call("rbind", dat)
names(dat) <- c("ccode1","ccode2","Ss","year")

tmp <- lapply(Ts.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  tmp[[w]]$year <- yr[w]
}
tmp <- do.call("rbind", tmp)
names(tmp) <- c("ccode1","ccode2","Ts","year")

dat <- merge(dat,tmp); rm(tmp)

## Add all the covariates to the data frame,
## dyadic covariates first. Lag these all
## by a year.
for (w in d.cons) {
  tmp <- lapply(get(w), function(x) as.data.frame(as.table(x)))
  for (z in 1:length(tmp)) {
    tmp[[z]]$year <- yr[z] + 1
  }
  tmp <- do.call("rbind", tmp)
  names(tmp) <- c("ccode1","ccode2",w,"year")
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}; rm(w)

## Now merge monadic covariates.
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  if ( w == "dem" ) {
    dat[[paste(w,".both",sep="")]] <- dat[[paste(w,".cc1", sep="")]] * dat[[paste(w,".cc2", sep="")]]
  }
  rm(tmp)
}

for (w in b.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}

## Drop "self dyads
dat <- dat[dat$ccode1 != dat$ccode2,]

####Make the monadic data for models of intrastate conflict
dat_mon <-mons.total[,c("ccode1","year","intra_intens")]
names(dat_mon) <- c("ccode1","year","intra_intens")
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",w)
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat_mon <- merge(dat_mon,tmp,all.x=TRUE)
  rm(tmp)
}


## Make unique identifiers
dat$dyadid <- (as.numeric(dat$ccode1) * 1000) + as.numeric(dat$ccode2)

## Now collect all the nonmissing observations
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
## Specify training set and validation set
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)



##### Estimate logit models
logit.threatn.train <- glm(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.threatn.train, dt[idx.out,], type="response", )
pred.use.threatn.08 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.threatn.08) <- c("observe","ccode1","ccode2","predict")

logit.basic.threatn.train <- glm(Ts ~ intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.basic.threatn.train, dt[idx.out,], type="response", )
pred.use.basic.threatn.08 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.basic.threatn.08) <- c("observe","ccode1","ccode2","predict")


## Now do the same for Intrastate conflict, constraining intra_intens to binary
## Specify training set and validation set
dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)

logit.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.intran.08 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.intran.08) <- c("observe","ccode1","predict")

logit.basic.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.basic.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.basic.intran.08 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.basic.intran.08) <- c("observe","ccode1","predict")


###############################
##2009
###############################

keep(dir, n.clus, mons.total, nets.total, trainyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, saom.pre.threat.net.07, saom.pre.intra_intens.07, pred.use.threatn.07, pred.use.intran.07, pred.use.basic.threatn.07, pred.use.basic.intran.07, saom.pre.threat.net.08, saom.pre.intra_intens.08, pred.use.threatn.08, pred.use.intran.08, pred.use.basic.threatn.08, pred.use.basic.intran.08, sure=TRUE)
setwd(dir)

testyr<-2009

## Specify data for out-of-sample prediction
yr <- seq((testyr-1), testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Must specify all covariates as constant (i.e., not changing over time)
dyad.effects <- d.cons
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}; rm(x)

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), coCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
  
)

eff <- getEffects(netdata)

###Specify the SAOM model

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect

eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add information from the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")


## Grab estimates from training sample
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.train.paper_04")

## Save rate and degree start values, since these must differ
## from training period.
rp1 <- eff[eff$effectName=="basic rate parameter support.net",]$initialValue
rp2 <- eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp3 <- eff[eff$effectName=="basic rate parameter threat.net",]$initialValue
rp4 <- eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp5 <- eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue
rp6 <- eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue

## Update thetas with estimates from training period
eff <- updateTheta(eff,output)
## Then replace rate and density parameters
eff[eff$effectName=="basic rate parameter support.net",]$initialValue <- rp1
eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp2
eff[eff$effectName=="basic rate parameter threat.net",]$initialValue <- rp3
eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp4
eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue <- rp5
eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue <- rp6

## No subroutines in stage 2, since we're not estimating parameters. Instead,
## just simulate the network a few thousand times.
model.SIM <- sienaModelCreate(useStdInits = F,
                              projname="Simulation-GOF-Predict-09",
                              nsub=0, n3=3000, maxlike=F, cond=F, modelType=c(support.net=1,threat.net=1))

## Call the estimation/simulation function. Be sure to return the dependent
## simulated matrices.
(n.clus <- detectCores() - 1)
SIM.out.predict <- siena07(model.SIM, data=netdata, effects=eff,
                           batch=TRUE, verbose=FALSE, useCluster=TRUE,
                           nbrNodes=n.clus, initC=TRUE, returnDeps=TRUE)
save(SIM.out.predict,file=paste("output.GOF.predict.paper_09"))


keep(mons.total, nets.total, dir, n.clus, id, N, yr, Ss.NET, Ts.NET, intra_intens, trainyr, testyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, saom.pre.threat.net.07, saom.pre.intra_intens.07, pred.use.threatn.07, pred.use.intran.07, pred.use.basic.threatn.07, pred.use.basic.intran.07, saom.pre.threat.net.08, saom.pre.intra_intens.08, pred.use.threatn.08, pred.use.intran.08, pred.use.basic.threatn.08, pred.use.basic.intran.08, sure=TRUE)

## Load SAOM predictions
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.predict.paper_09")



## Collect SAOM predictions for threat and intrastate conflict
sims <- SIM.out.predict$sims

#Start for threat ties
## Grab predicted SAOM networks for last time period
base <- matrix(0,N,N,dimnames=list(id,id))
for ( i in 1:length(sims) ) {
  adj <- matrix(0,N,N)
  ## The final year is the one we're predicting
  edges <- sims[[i]][[1]][["threat.net"]][[(length(yr)-1)]]
  adj[edges[, 1:2]] <- edges[,3]
  base <- base + adj
  rm(adj,edges)
}
## Weight total predicted ties by number of simulated networks
base <- base/length(sims); diag(base) <- NA

## Get the * observed* networks for last time period
obs <- Ts.NET[[length(yr)]]
diag(obs) <- NA

## Combine predicted and observed vectors
saom.pre.threat.net.09 <- na.omit(cbind(as.data.frame(as.table(base)),as.data.frame(as.table(obs))[,3]))
names(saom.pre.threat.net.09) <- c("ccode1","ccode2","predict","observe")    

#Now for Intrastate violence -- specifically looking at all armed conflict
base<-as.vector(matrix(0,nrow=N,dimnames=list(id)))
for ( i in 1:length(sims) ) {
  ## The final year is the one we're predicting
  ac <- sims[[i]][[1]][["intra_intens.net"]][[(length(yr)-1)]]
  ac[ac>=1]<-1
  base <- base + ac
  rm(ac)
}

## Weight total predicted ties by number of simulated networks
base <- base/length(sims)

## Get the observed civil wars from last time period
intra.obs<-intra_intens[,length(yr)]

##combine predicted and observed
op <- na.omit(cbind(base,intra.obs,id))
saom.pre.intra_intens.09<-as.data.frame(op)
names(saom.pre.intra_intens.09) <- c("predict","observe","ccode1")
saom.pre.intra_intens.09$observe[saom.pre.intra_intens.09$observe>=1]<-1

rm(base,op,sims)

#################################################

## Now estimate a logit model, predict out of sample, and compare
## the logit to the SAOM predictions

## Reorganize the dataset for the logit model. Use the period up until the test year.
## We're going to grab all the same data as for the SAOMs, then
## restructure it into a data frame.
yr <- seq(1980, testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with the DVs
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))

#########################################################
## Then do all the covariates for both networks
#########################################################


## Names of dyadic and monadic covariates, plus behavior variable
#Create new name for the intra_intens that will be a covariate
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]
mons.total$intra_intensc<-mons.total$intra_intens
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year") ]


## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################################
## Now convert all the data matrices to a data frame - one dyadic and one monadic
####################################################

dat <- lapply(Ss.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  dat[[w]]$year <- yr[w]
}
dat <- do.call("rbind", dat)
names(dat) <- c("ccode1","ccode2","Ss","year")

tmp <- lapply(Ts.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  tmp[[w]]$year <- yr[w]
}
tmp <- do.call("rbind", tmp)
names(tmp) <- c("ccode1","ccode2","Ts","year")

dat <- merge(dat,tmp); rm(tmp)

## Add all the covariates to the data frame,
## dyadic covariates first. Lag these all
## by a year.
for (w in d.cons) {
  tmp <- lapply(get(w), function(x) as.data.frame(as.table(x)))
  for (z in 1:length(tmp)) {
    tmp[[z]]$year <- yr[z] + 1
  }
  tmp <- do.call("rbind", tmp)
  names(tmp) <- c("ccode1","ccode2",w,"year")
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}; rm(w)

## Now merge monadic covariates.
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  if ( w == "dem" ) {
    dat[[paste(w,".both",sep="")]] <- dat[[paste(w,".cc1", sep="")]] * dat[[paste(w,".cc2", sep="")]]
  }
  rm(tmp)
}

for (w in b.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}

## Drop "self dyads
dat <- dat[dat$ccode1 != dat$ccode2,]

####Make the monadic data for models of intrastate conflict
dat_mon <-mons.total[,c("ccode1","year","intra_intens")]
names(dat_mon) <- c("ccode1","year","intra_intens")
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",w)
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat_mon <- merge(dat_mon,tmp,all.x=TRUE)
  rm(tmp)
}


## Make unique identifiers
dat$dyadid <- (as.numeric(dat$ccode1) * 1000) + as.numeric(dat$ccode2)

## Now collect all the nonmissing observations
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
## Specify training set and validation set
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)



##### Estimate logit models
logit.threatn.train <- glm(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.threatn.train, dt[idx.out,], type="response", )
pred.use.threatn.09 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.threatn.09) <- c("observe","ccode1","ccode2","predict")

logit.basic.threatn.train <- glm(Ts ~ intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.basic.threatn.train, dt[idx.out,], type="response", )
pred.use.basic.threatn.09 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.basic.threatn.09) <- c("observe","ccode1","ccode2","predict")


## Now do the same for Intrastate conflict, constraining intra_intens to binary
## Specify training set and validation set
dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)

logit.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.intran.09 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.intran.09) <- c("observe","ccode1","predict")

logit.basic.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.basic.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.basic.intran.09 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.basic.intran.09) <- c("observe","ccode1","predict")




###############################
##2010
###############################

keep(dir, n.clus, mons.total, nets.total, trainyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, saom.pre.threat.net.07, saom.pre.intra_intens.07, pred.use.threatn.07, pred.use.intran.07, pred.use.basic.threatn.07, pred.use.basic.intran.07, saom.pre.threat.net.08, saom.pre.intra_intens.08, pred.use.threatn.08, pred.use.intran.08, pred.use.basic.threatn.08, pred.use.basic.intran.08,  saom.pre.threat.net.09, saom.pre.intra_intens.09, pred.use.threatn.09, pred.use.intran.09, pred.use.basic.threatn.09, pred.use.basic.intran.09, sure=TRUE)
setwd(dir)

testyr<-2010

## Specify data for out-of-sample prediction
yr <- seq((testyr-1), testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, i.e., for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Must specify all covariates as constant (i.e., not changing over time)
dyad.effects <- d.cons
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}; rm(x)

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), coCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
  
)

eff <- getEffects(netdata)

###Specify the SAOM model

## Support network equation first
d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff <- includeEffects(eff, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff <- includeEffects(eff, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff <- includeEffects(eff, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff <- includeEffects(eff, recip, name="support.net", include=T, type="eval")

## Now the threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Include the cross-products
eff <- includeEffects(eff, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff <- includeEffects(eff, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff <- includeEffects(eff, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")

## Add friend of my enemy effect

eff <- includeEffects(eff, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add information from the support groups
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff <- includeEffects(eff, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff <- includeInteraction(eff, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")


## Grab estimates from training sample
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.train.paper_04")

## Save rate and degree start values, since these must differ
## from training period.
rp1 <- eff[eff$effectName=="basic rate parameter support.net",]$initialValue
rp2 <- eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp3 <- eff[eff$effectName=="basic rate parameter threat.net",]$initialValue
rp4 <- eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue
rp5 <- eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue
rp6 <- eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue

## Update thetas with estimates from training period
eff <- updateTheta(eff,output)
## Then replace rate and density parameters
eff[eff$effectName=="basic rate parameter support.net",]$initialValue <- rp1
eff[eff$effectName=="support.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp2
eff[eff$effectName=="basic rate parameter threat.net",]$initialValue <- rp3
eff[eff$effectName=="threat.net: outdegree (density)" & eff$include==TRUE,]$initialValue <- rp4
eff[eff$effectName=="rate intra_intens.net period 1",]$initialValue <- rp5
eff[eff$effectName=="intra_intens.net linear shape" & eff$include==TRUE,]$initialValue <- rp6

## No subroutines in stage 2, since we're not estimating parameters. Instead,
## just simulate the network a few thousand times.
model.SIM <- sienaModelCreate(useStdInits = F,
                              projname="Simulation-GOF-Predict-10",
                              nsub=0, n3=3000, maxlike=F, cond=F, modelType=c(support.net=1,threat.net=1))

## Call the estimation/simulation function. Be sure to return the dependent
## simulated matrices.
(n.clus <- detectCores() - 1)
SIM.out.predict <- siena07(model.SIM, data=netdata, effects=eff,
                           batch=TRUE, verbose=FALSE, useCluster=TRUE,
                           nbrNodes=n.clus, initC=TRUE, returnDeps=TRUE)
save(SIM.out.predict,file=paste("output.GOF.predict.paper_10"))


keep(mons.total, nets.total, dir, n.clus, id, N, yr, Ss.NET, Ts.NET, intra_intens, trainyr, testyr, saom.pre.threat.net.05, saom.pre.intra_intens.05, pred.use.threatn.05, pred.use.intran.05, pred.use.basic.threatn.05, pred.use.basic.intran.05, saom.pre.threat.net.06, saom.pre.intra_intens.06, pred.use.threatn.06, pred.use.intran.06, pred.use.basic.threatn.06, pred.use.basic.intran.06, saom.pre.threat.net.07, saom.pre.intra_intens.07, pred.use.threatn.07, pred.use.intran.07, pred.use.basic.threatn.07, pred.use.basic.intran.07, saom.pre.threat.net.08, saom.pre.intra_intens.08, pred.use.threatn.08, pred.use.intran.08, pred.use.basic.threatn.08, pred.use.basic.intran.08, saom.pre.threat.net.09, saom.pre.intra_intens.09, pred.use.threatn.09, pred.use.intran.09, pred.use.basic.threatn.09, pred.use.basic.intran.09, sure=TRUE)

## Load SAOM predictions
setwd(paste(dir, "/Output", sep=""))
load("output.GOF.predict.paper_10")



## Collect SAOM predictions for threat and intrastate conflict
sims <- SIM.out.predict$sims

#Start for threat ties
## Grab predicted SAOM networks for last time period
base <- matrix(0,N,N,dimnames=list(id,id))
for ( i in 1:length(sims) ) {
  adj <- matrix(0,N,N)
  ## The final year is the one we're predicting
  edges <- sims[[i]][[1]][["threat.net"]][[(length(yr)-1)]]
  adj[edges[, 1:2]] <- edges[,3]
  base <- base + adj
  rm(adj,edges)
}
## Weight total predicted ties by number of simulated networks
base <- base/length(sims); diag(base) <- NA

## Get the * observed* networks for last time period
obs <- Ts.NET[[length(yr)]]
diag(obs) <- NA

## Combine predicted and observed vectors
saom.pre.threat.net.10 <- na.omit(cbind(as.data.frame(as.table(base)),as.data.frame(as.table(obs))[,3]))
names(saom.pre.threat.net.10) <- c("ccode1","ccode2","predict","observe")    

#Now for Intrastate violence -- specifically looking at all armed conflict
base<-as.vector(matrix(0,nrow=N,dimnames=list(id)))
for ( i in 1:length(sims) ) {
  ## The final year is the one we're predicting
  ac <- sims[[i]][[1]][["intra_intens.net"]][[(length(yr)-1)]]
  ac[ac>=1]<-1
  base <- base + ac
  rm(ac)
}

## Weight total predicted ties by number of simulated networks
base <- base/length(sims)

## Get the observed civil wars from last time period
intra.obs<-intra_intens[,length(yr)]

##combine predicted and observed
op <- na.omit(cbind(base,intra.obs,id))
saom.pre.intra_intens.10<-as.data.frame(op)
names(saom.pre.intra_intens.10) <- c("predict","observe","ccode1")
saom.pre.intra_intens.10$observe[saom.pre.intra_intens.10$observe>=1]<-1

rm(base,op,sims)

#################################################

## Now estimate a logit model, predict out of sample, and compare
## the logit to the SAOM predictions

## Reorganize the dataset for the logit model. Use the period up until the test year.
## We're going to grab all the same data as for the SAOMs, then
## restructure it into a data frame.
yr <- seq(1980, testyr)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

#########################################################
## Start with the DVs
#########################################################


##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))

#########################################################
## Then do all the covariates for both networks
#########################################################


## Names of dyadic and monadic covariates, plus behavior variable
#Create new name for the intra_intens that will be a covariate
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]
mons.total$intra_intensc<-mons.total$intra_intens
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year") ]


## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not subtracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################################
## Now convert all the data matrices to a data frame - one dyadic and one monadic
####################################################

dat <- lapply(Ss.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  dat[[w]]$year <- yr[w]
}
dat <- do.call("rbind", dat)
names(dat) <- c("ccode1","ccode2","Ss","year")

tmp <- lapply(Ts.NET, function(x) as.data.frame(as.table(x)))
for (w in 1:length(yr)) {
  tmp[[w]]$year <- yr[w]
}
tmp <- do.call("rbind", tmp)
names(tmp) <- c("ccode1","ccode2","Ts","year")

dat <- merge(dat,tmp); rm(tmp)

## Add all the covariates to the data frame,
## dyadic covariates first. Lag these all
## by a year.
for (w in d.cons) {
  tmp <- lapply(get(w), function(x) as.data.frame(as.table(x)))
  for (z in 1:length(tmp)) {
    tmp[[z]]$year <- yr[z] + 1
  }
  tmp <- do.call("rbind", tmp)
  names(tmp) <- c("ccode1","ccode2",w,"year")
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}; rm(w)

## Now merge monadic covariates.
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  if ( w == "dem" ) {
    dat[[paste(w,".both",sep="")]] <- dat[[paste(w,".cc1", sep="")]] * dat[[paste(w,".cc2", sep="")]]
  }
  rm(tmp)
}

for (w in b.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",paste(w,".cc1", sep=""))
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat <- merge(dat,tmp,all.x=TRUE)
  names(tmp) <- c("ccode2","year",paste(w,".cc2", sep=""))
  dat <- merge(dat,tmp,all.x=TRUE)
  rm(tmp)
}

## Drop "self dyads
dat <- dat[dat$ccode1 != dat$ccode2,]

####Make the monadic data for models of intrastate conflict
dat_mon <-mons.total[,c("ccode1","year","intra_intens")]
names(dat_mon) <- c("ccode1","year","intra_intens")
for (w in m.cons) {
  tmp <- as.data.frame(as.table(get(w)))
  names(tmp) <- c("ccode1","year",w)
  tmp$year <- as.numeric(as.character(tmp$year)) + 1
  dat_mon <- merge(dat_mon,tmp,all.x=TRUE)
  rm(tmp)
}


## Make unique identifiers
dat$dyadid <- (as.numeric(dat$ccode1) * 1000) + as.numeric(dat$ccode2)

## Now collect all the nonmissing observations
dt <- na.omit(dat[,c("Ss", "Ts", "threat_density.cc1", "threat_recip", "support_mi.cc1", "support_mi.cc2", "support_ipcloseness.cc1", "support_ipcloseness.cc2", "intra_intens.cc1", "intra_intens.cc2", "IGOcommon", "idealpointdistance", "contiguous", "dem.cc1", "dem.cc2", "dem.both", "mi1.cc1", "mi1.cc2", "energypc.cc1", "energypc.cc2", "mmp_change.cc1", "mmp_change.cc2", "year","dyadid","ccode1","ccode2")])
dt <- dt[order(dt$dyadid,dt$year),]
## Specify training set and validation set
idx.in <- which(dt$year <= trainyr)
idx.out <- which(dt$year == testyr)



##### Estimate logit models
logit.threatn.train <- glm(Ts ~ Ss + intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip + support_mi.cc1 * support_ipcloseness.cc1 + support_mi.cc2 * support_ipcloseness.cc2, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.threatn.train, dt[idx.out,], type="response", )
pred.use.threatn.10 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.threatn.10) <- c("observe","ccode1","ccode2","predict")

logit.basic.threatn.train <- glm(Ts ~ intra_intens.cc1 + intra_intens.cc2 + IGOcommon + idealpointdistance + contiguous + dem.cc1 + dem.cc2 + dem.both + mi1.cc1 + mi1.cc2 + energypc.cc1 + energypc.cc2 + mmp_change.cc1 + mmp_change.cc2 + threat_recip, data=dt, family=binomial(link="logit"), subset=idx.in)

## Now use the above estimated model to make predictions out of sample
tst <- predict(logit.basic.threatn.train, dt[idx.out,], type="response", )
pred.use.basic.threatn.10 <- cbind( dt[dt$year==testyr, c("Ts","ccode1","ccode2")], tst)
names(pred.use.basic.threatn.10) <- c("observe","ccode1","ccode2","predict")


## Now do the same for Intrastate conflict, constraining intra_intens to binary
## Specify training set and validation set
dt_m <- na.omit(dat_mon[,c("intra_intens", "intra_intensc", "dem", "mi1","energypc", "mmp_change", "support_mi", "threat_mi", "support_ipcloseness", "year","ccode1")])
dt_m <- dt_m[order(dt_m$ccode1,dt_m$year),]
dt_m$intra_intens[dt_m$intra_intens>=1]<-1
idx.in <- which(dt_m$year <= trainyr)
idx.out <- which(dt_m$year == testyr)

logit.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change + threat_mi + support_mi * support_ipcloseness, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.intran.10 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.intran.10) <- c("observe","ccode1","predict")

logit.basic.intran.train <- glm(intra_intens ~ intra_intensc + I(intra_intensc^2) + dem + mi1 + energypc + mmp_change, data=dt_m, family=binomial(link="logit"), subset=idx.in)

tst <- predict(logit.basic.intran.train, dt_m[idx.out,], type=c("response"))
pred.use.basic.intran.10 <- cbind(dt_m[dt_m$year==testyr,c("intra_intens","ccode1")], tst)
names(pred.use.basic.intran.10) <- c("observe","ccode1","predict")



######################################################################
###Pull together the prediction data frames
saom.pre.threat.net.05<-data.frame(saom.pre.threat.net.05, year=2005)
saom.pre.threat.net.06<-data.frame(saom.pre.threat.net.06, year=2006)
saom.pre.threat.net.07<-data.frame(saom.pre.threat.net.07, year=2007)
saom.pre.threat.net.08<-data.frame(saom.pre.threat.net.08, year=2008)
saom.pre.threat.net.09<-data.frame(saom.pre.threat.net.09, year=2009)
saom.pre.threat.net.10<-data.frame(saom.pre.threat.net.10, year=2010)
saom.pre.intra_intens.05<-data.frame(saom.pre.intra_intens.05, year=2005)
saom.pre.intra_intens.06<-data.frame(saom.pre.intra_intens.06, year=2006)
saom.pre.intra_intens.07<-data.frame(saom.pre.intra_intens.07, year=2007)
saom.pre.intra_intens.08<-data.frame(saom.pre.intra_intens.08, year=2008)
saom.pre.intra_intens.09<-data.frame(saom.pre.intra_intens.09, year=2009)
saom.pre.intra_intens.10<-data.frame(saom.pre.intra_intens.10, year=2010)
pred.use.threatn.05<-data.frame(pred.use.threatn.05, year=2005)
pred.use.threatn.06<-data.frame(pred.use.threatn.06, year=2006)
pred.use.threatn.07<-data.frame(pred.use.threatn.07, year=2007)
pred.use.threatn.08<-data.frame(pred.use.threatn.08, year=2008)
pred.use.threatn.09<-data.frame(pred.use.threatn.09, year=2009)
pred.use.threatn.10<-data.frame(pred.use.threatn.10, year=2010)
pred.use.intran.05<-data.frame(pred.use.intran.05, year=2005)
pred.use.intran.06<-data.frame(pred.use.intran.06, year=2006)
pred.use.intran.07<-data.frame(pred.use.intran.07, year=2007)
pred.use.intran.08<-data.frame(pred.use.intran.08, year=2008)
pred.use.intran.09<-data.frame(pred.use.intran.09, year=2009)
pred.use.intran.10<-data.frame(pred.use.intran.10, year=2010)
pred.use.basic.threatn.05<-data.frame(pred.use.basic.threatn.05, year=2005)
pred.use.basic.threatn.06<-data.frame(pred.use.basic.threatn.06, year=2006)
pred.use.basic.threatn.07<-data.frame(pred.use.basic.threatn.07, year=2007)
pred.use.basic.threatn.08<-data.frame(pred.use.basic.threatn.08, year=2008)
pred.use.basic.threatn.09<-data.frame(pred.use.basic.threatn.09, year=2009)
pred.use.basic.threatn.10<-data.frame(pred.use.basic.threatn.10, year=2010)
pred.use.basic.intran.05<-data.frame(pred.use.basic.intran.05, year=2005)
pred.use.basic.intran.06<-data.frame(pred.use.basic.intran.06, year=2006)
pred.use.basic.intran.07<-data.frame(pred.use.basic.intran.07, year=2007)
pred.use.basic.intran.08<-data.frame(pred.use.basic.intran.08, year=2008)
pred.use.basic.intran.09<-data.frame(pred.use.basic.intran.09, year=2009)
pred.use.basic.intran.10<-data.frame(pred.use.basic.intran.10, year=2010)
##################################################################################################################################
#Do figures
##################################################################################################################################

## Create the figures comparing SAOM to Logit; Appendix Figure 12
library(ggplot2)
library(grid)
library(gridExtra)

## Append the prediciton files
setwd(paste(dir, "/Output", sep=""))

saom.pre.threat.net<-rbind(saom.pre.threat.net.05, saom.pre.threat.net.06, saom.pre.threat.net.07, saom.pre.threat.net.08, saom.pre.threat.net.09, saom.pre.threat.net.10)
write.csv(saom.pre.threat.net, file = "saom.pre.threat.net_paper.csv", row.names=FALSE)
saom.pre.intra_intens<-rbind(saom.pre.intra_intens.05, saom.pre.intra_intens.06, saom.pre.intra_intens.07, saom.pre.intra_intens.08, saom.pre.intra_intens.09, saom.pre.intra_intens.10)
write.csv(saom.pre.intra_intens, file = "saom.pre.intra_intens_paper.csv", row.names=FALSE)
pred.use.threatn<-rbind(pred.use.threatn.05, pred.use.threatn.06, pred.use.threatn.07, pred.use.threatn.08, pred.use.threatn.09, pred.use.threatn.10)
write.csv(pred.use.threatn, file = "pred.use.threatn_paper.csv", row.names=FALSE)
pred.use.intran<-rbind(pred.use.intran.05, pred.use.intran.06, pred.use.intran.07, pred.use.intran.08, pred.use.intran.09, pred.use.intran.10)
write.csv(pred.use.intran, file = "pred.use.intran_paper.csv", row.names=FALSE)
pred.use.basic.threatn<-rbind(pred.use.basic.threatn.05, pred.use.basic.threatn.06, pred.use.basic.threatn.07, pred.use.basic.threatn.08, pred.use.basic.threatn.09, pred.use.basic.threatn.10)
write.csv(pred.use.basic.threatn, file = "pred.use.basic.threatn_paper.csv", row.names=FALSE)
pred.use.basic.intran<-rbind(pred.use.basic.intran.05, pred.use.basic.intran.06, pred.use.basic.intran.07, pred.use.basic.intran.08, pred.use.basic.intran.09, pred.use.basic.intran.10)
write.csv(pred.use.basic.intran, file = "pred.use.basic.intran_paper.csv", row.names=FALSE)

keep(dir, n.clus, sure=TRUE)

setwd(paste(dir, "/Output", sep=""))

saom.pre.threat.net <- read.csv("saom.pre.threat.net_paper.csv", header=TRUE, row.names=NULL)
saom.pre.intra_intens <- read.csv("saom.pre.intra_intens_paper.csv", header=TRUE, row.names=NULL)
pred.use.threatn <- read.csv("pred.use.threatn_paper.csv", header=TRUE, row.names=NULL)
pred.use.intran <- read.csv("pred.use.intran_paper.csv", header=TRUE, row.names=NULL)
#pred.use.basic.threatn <- read.csv("pred.use.basic.threatn_paper.csv", header=TRUE, row.names=NULL)
#pred.use.basic.intran <- read.csv("pred.use.basic.intran_paper.csv", header=TRUE, row.names=NULL)

## Who are the most active threat countries in 2005-2010?
threats.ct <- aggregate(saom.pre.threat.net$observe, by=list(Category=saom.pre.threat.net$ccode1), FUN=sum)
threats.focus <- threats.ct$Category[threats.ct$x>=18]

## Who has intrastate conflict in 2005-2010?
intra.ct <- aggregate(saom.pre.intra_intens$observe, by=list(Category=saom.pre.intra_intens$ccode1), FUN=sum)
intra.focus <- intra.ct$Category[intra.ct$x==6]
rm(threats.ct, intra.ct)

## Find cumulative predicted vs. observed threats for focus states with threats, first for SAOMs
use.threat.saom <- na.omit(saom.pre.threat.net[saom.pre.threat.net$ccode1 %in% threats.focus,])
threat.pred_obs.saom <- data.frame(
    ccode1 = aggregate(use.threat.saom$predict, by=list(Category=use.threat.saom$ccode1), FUN=sum)$Category,
    pred_obs = aggregate(use.threat.saom$predict,
        by=list(Category=use.threat.saom$ccode1)
        , FUN=sum)$x / (
              aggregate(use.threat.saom$observe,
                        by=list(Category=use.threat.saom$ccode1),
                        FUN=sum)$x 
              )
    )
rm(use.threat.saom)

## Then for threats logit model

use.threatn <- na.omit(pred.use.threatn[pred.use.threatn$ccode1 %in% threats.focus,])
threatn.pred_obs.logit <- data.frame(
  ccode1 = aggregate(use.threatn$predict, by=list(Category=use.threatn$ccode1), FUN=sum)$Category,
  pred_obs = aggregate(use.threatn$predict,
                       by=list(Category=use.threatn$ccode1)
                       , FUN=sum)$x / (
                         aggregate(use.threatn$observe,
                                   by=list(Category=use.threatn$ccode1),
                                   FUN=sum)$x 
                       )
)

rm(use.threatn)


## Do the same for intrastate violence
use.intra.saom <- na.omit(saom.pre.intra_intens[saom.pre.intra_intens$ccode1 %in% intra.focus,])
intra.pred_obs.saom<-data.frame(
  ccode1 = aggregate(use.intra.saom$predict, by=list(Category=use.intra.saom$ccode1), FUN=sum)$Category,
  pred_obs = aggregate(use.intra.saom$predict,
                       by=list(Category=use.intra.saom$ccode1)
                       , FUN=sum)$x / (
                         aggregate(use.intra.saom$observe,
                                   by=list(Category=use.intra.saom$ccode1),
                                   FUN=sum)$x 
                       )
)
rm(use.intra.saom)

## Now logit

use.intran <- na.omit(pred.use.intran[pred.use.intran$ccode1 %in% intra.focus,])
intran.pred_obs.logit<-data.frame(
  ccode1 = aggregate(use.intran$predict, by=list(Category=use.intran$ccode1), FUN=sum)$Category,
  pred_obs = aggregate(use.intran$predict,
                       by=list(Category=use.intran$ccode1)
                       , FUN=sum)$x / (
                         aggregate(use.intran$observe,
                                   by=list(Category=use.intran$ccode1),
                                   FUN=sum)$x 
                       )
)

rm(use.intran)


## Plot the predictions in a barplot
df.pred_obs.threat <- rbind(threatn.pred_obs.logit,threat.pred_obs.saom)
df.pred_obs.threat[df.pred_obs.threat$pred_obs<0.01,"pred_obs"] <- 0.01 # Set zeros to negligibly small value for plotting
df.pred_obs.threat$model <- c(rep("Logit",12),rep("SAOM",12))
df.pred_obs.threat$cow <- countrycode(df.pred_obs.threat$ccode1, "cown", "cowc", warn = FALSE)
ord <- df.pred_obs.threat[df.pred_obs.threat$model=="SAOM",][order(df.pred_obs.threat[df.pred_obs.threat$model=="SAOM",]$pred_obs),]$cow
df.pred_obs.threat$cow <- factor(df.pred_obs.threat$cow, levels=ord)

gg.pred_obs.threat <- ggplot(df.pred_obs.threat, aes(x=cow, y=pred_obs, group=model, fill=model)) +
    geom_bar(stat="identity", position="dodge", width=.65) +
    scale_fill_manual("", values=c("#0072B2", "#8519d7"), labels=c("Logit", "SAOM")) +
theme(
    legend.position="bottom",
    legend.margin=margin(0,0,0,0,"mm"),
    axis.text.x=element_text(color="black", angle=45, hjust=1),
    plot.title = element_text(hjust = 0.5)
    ) +
    ylab("Ratio of Predicted vs. Observed") +
    xlab("") +
    ggtitle("Interstate threats")

## Now for intrastate war
df.pred_obs.intra <- rbind(intran.pred_obs.logit, intra.pred_obs.saom)
df.pred_obs.intra[df.pred_obs.intra$pred_obs<0.01,"pred_obs"] <- 0.01 # Set zeros to negligibly small value for plotting
df.pred_obs.intra$model <- c(rep("Logit",17),rep("SAOM",17))
df.pred_obs.intra$cow <- countrycode(df.pred_obs.intra$ccode1, "cown", "cowc", warn = FALSE)
ord <- df.pred_obs.intra[df.pred_obs.intra$model=="SAOM",][order(df.pred_obs.intra[df.pred_obs.intra$model=="SAOM",]$pred_obs),]$cow
df.pred_obs.intra$cow <- factor(df.pred_obs.intra$cow, levels=ord)

gg.pred_obs.intra <- ggplot(df.pred_obs.intra, aes(x=cow, y=pred_obs, group=model, fill=model)) +
    geom_bar(stat="identity", position="dodge", width=.65) +
    scale_fill_manual("", values=c("#0072B2", "#8519d7"), labels=c("Logit", "SAOM")) +
theme(
    legend.position="bottom",
    legend.margin=margin(0,0,0,0,"mm"),
    axis.text.x=element_text(color="black", angle=45, hjust=1),
    plot.title = element_text(hjust = 0.5)
    ) +
    ylab("Ratio of Predicted vs. Observed") +
    xlab("") +
    ggtitle("Intrastate Armed Conflict")



setwd(paste(dir, "/Figures", sep=""))

pdf("pred_obs-paper_threat.pdf", width=5, height=3.5, paper="special")
par(mar=c(0,0,0,0))
gg.pred_obs.threat
dev.off()

pdf("pred_obs-paper_intra.pdf", width=5, height=3.5, paper="special")
par(mar=c(0,0,0,0))
gg.pred_obs.intra
dev.off()



#########Generate ROC and PR curves for the sets of cases with threats or intrastate armed conflict
###Appendix Figure 13

#First subset to cases with any threats or intrastate armed conflict during the entire data
setwd(dir)
nets.total <- read.csv("relation_network.csv", header=TRUE, row.names=1)
nets.total$dyadid <- (as.numeric(nets.total$ccode1) * 1000) + as.numeric(nets.total$ccode2)
threats.ct <- aggregate(nets.total$threat_bin_3, by=list(Category=nets.total$dyadid), FUN=sum)
threats.any <- threats.ct$Category[threats.ct$x>=1]

mons.total <- read.csv("relation_nodal.csv", header=TRUE, row.names=1)
intra.ct <- aggregate(mons.total$intra_intens, by=list(Category=mons.total$ccode1), FUN=sum)
intra.any <- intra.ct$Category[intra.ct$x>=1]

saom.pre.threat.net$dyadid<-(as.numeric(saom.pre.threat.net$ccode1) * 1000) + as.numeric(saom.pre.threat.net$ccode2)
use.threat.saom <- na.omit(saom.pre.threat.net[saom.pre.threat.net$dyadid %in% threats.any,])
pred.use.threatn$dyadid<-(as.numeric(pred.use.threatn$ccode1) * 1000) + as.numeric(pred.use.threatn$ccode2)
use.threatn.logit <- na.omit(pred.use.threatn[pred.use.threatn$dyadid %in% threats.any,])

use.intra.saom <- na.omit(saom.pre.intra_intens[saom.pre.intra_intens$ccode1 %in% intra.any,])
use.intran.logit <- na.omit(pred.use.intran[pred.use.intran$ccode1 %in% intra.any,])

#Generate prediction estimates
pred.threat.saom <- prediction(use.threat.saom$predict, use.threat.saom$observe)
pred.threatn.logit <- prediction(use.threatn.logit$predict, use.threatn.logit$observe)
pred.intra.saom <- prediction(use.intra.saom$predict, use.intra.saom$observe)
pred.intran.logit <- prediction(use.intran.logit$predict, use.intran.logit$observe)

roc.threat.saom <- performance( pred.threat.saom, "tpr", "fpr" )
prc.threat.saom <- performance( pred.threat.saom, "prec", "rec" )
auc.threat.saom <- performance( pred.threat.saom, "auc")
aucpr.threat.saom <- performance( pred.threat.saom, "aucpr")

roc.threatn.logit <- performance( pred.threatn.logit, "tpr", "fpr" )
prc.threatn.logit <- performance( pred.threatn.logit, "prec", "rec" )
auc.threatn.logit <- performance( pred.threatn.logit, "auc")
aucpr.threatn.logit <- performance( pred.threatn.logit, "aucpr")

roc.intra.saom <- performance( pred.intra.saom, "tpr", "fpr" )
prc.intra.saom <- performance( pred.intra.saom, "prec", "rec" )
auc.intra.saom <- performance( pred.intra.saom, "auc")
aucpr.intra.saom <- performance( pred.intra.saom, "aucpr")

roc.intran.logit <- performance( pred.intran.logit, "tpr", "fpr" )
prc.intran.logit <- performance( pred.intran.logit, "prec", "rec" )
auc.intran.logit <- performance( pred.intran.logit, "auc")
aucpr.intran.logit <- performance( pred.intran.logit, "aucpr")


setwd(paste(dir, "/Figures", sep=""))

pdf("roc-paper_threat.pdf", width=5, height=3.5, paper="special")
plot( roc.threat.saom, col="#8519d7", lwd=3, main="ROC Curves, Interstate Threats" )
plot( roc.threatn.logit, add = TRUE, col = "#0072B2", lwd=3)
dev.off()

pdf("prc-paper_threat.pdf", width=5, height=3.5, paper="special")
plot( prc.threat.saom, col="#8519d7", lwd=3, main="Precision-Recall Curves, Interstate Threats" )
plot( prc.threatn.logit, add = TRUE, col = "#0072B2", lwd=3)
dev.off()

pdf("roc-paper_intra.pdf", width=5, height=3.5, paper="special")
plot( roc.intran.logit, col = "#0072B2", lwd=3, main="ROC Curves, Intrastate Armed Conflict")
plot( roc.intra.saom, add = TRUE, col="#8519d7", lwd=3 )
dev.off()

pdf("prc-paper_intra.pdf", width=5, height=3.5, paper="special")
plot( prc.intran.logit, col = "#0072B2", lwd=3, main="Precision-Recall Curves, Intrastate Armed Conflict")
plot( prc.intra.saom, add = TRUE, col="#8519d7", lwd=3 )
dev.off()



